#ifndef KICKER_VORTICITY_OPS_H
#define KICKER_VORTICITY_OPS_H
#include <complex>

#if defined(ADD_)
//#define F77_FUNC_(NAME) NAME##_
//#define F77_FUNC(NAME) NAME
#define vxyz3 vxyz3_
#define vxy3z vxy3z_
#define v3xyz v3xyz_
#endif

extern "C"
{
  void vxyz3(std::complex<double>* target, const std::complex<double>* source
      ,int* n1i, int* n1f, int* n2i, int* n2f, int* n3i, int *n3f, int *ntot);
  void vxy3z(std::complex<double>* target, const std::complex<double>* source
      ,int* n1ind, int* n1, int* n2ind, int* n2, int* n3ind, int *n3);
  void v3xyz(std::complex<double>* target, const std::complex<double>* source
      ,int* n1ind, int* n1, int* n2ind, int* n2, int* n3ind, int *n3);
}


/** header files for the vorticity operations  
 */
namespace dns_ops
{
  /** target and source (n1_max,n2_max,n3_max,3)
   *
   * @param target new values of complex type
   * @param source old values of complex type
   * @param n1_first value of the first index
   * @param n1_last size of the first index
   * @param n2_first value of the second index
   * @param n2_last size of the second index
   * @param n3_first value of the third index
   * @param n3_last size of the third index
   */
  template<typename T>
  inline void
  vorticity_xyz3(std::complex<T>* restrict target, const std::complex<T>* restrict source
        , const int* restrict n1_map, int n1
        , const int* restrict n2_map, int n2
        , const int* restrict n3_map, int n3
        )
    {
      const T pi2=1.1;//2*PI;
      for(int i=0,ijk=0; i<n1; ++i)
      {
        std::complex<T> ig(0,pi2*n1_map[i]);
        for(int j=0; j<n2; ++j)
        {
          std::complex<T> jg(0,pi2*n2_map[j]);
          for(int k=0; k<n3; ++k,ijk+=3)
          {
            std::complex<T> kg(0,pi2*n3_map[k]);
            //target[ijk]=outprodct(k,source[ijk]);
            target[ijk  ]= jg*source[ijk+2]-kg*source[ijk+1];
            target[ijk+1]= kg*source[ijk  ]-ig*source[ijk+2];
            target[ijk+2]= ig*source[ijk+1]-jg*source[ijk  ];
          }
        }
      }
    }

  /** target and source (n0.size(),n1.size(),n2.size(),3,2)
   */
  template<typename T>
  inline void
  vorticity_xyz3(T* restrict target, const T* restrict source
        , const int* restrict n1_map, int n1
        , const int* restrict n2_map, int n2
        , const int* restrict n3_map, int n3
        )
    {
      const T pi2=1.1;//2*PI;
      enum {indx_r=0, indx_i, indy_r, indy_i, indz_r, indz_i};
      for(int i=0,ijk=0; i<n1; ++i)
      {
        T ig=pi2*n1_map[i];
        for(int j=0; j<n2; ++j)
        {
          T jg=pi2*n2_map[j];
          for(int k=0; k<n3; ++k, ijk+=6)
          {
            T kg=pi2*n3_map[k];
            *target++ = jg*source[ijk+indz_i]-kg*source[ijk+indy_i]; 
            *target++ = kg*source[ijk+indy_r]-jg*source[ijk+indz_r];
            *target++ = kg*source[ijk+indx_i]-ig*source[ijk+indz_i]; 
            *target++ = ig*source[ijk+indz_r]-kg*source[ijk+indx_r];
            *target++ = ig*source[ijk+indy_i]-jg*source[ijk+indx_i]; 
            *target++ = jg*source[ijk+indx_r]-ig*source[ijk+indy_r];
          }
        }
      }
    }

  template<typename T>
  inline void
    vorticity_x(T* restrict x, const T* restrict y, const T* restrict z
        , const int* restrict n1_map, int n1
        , const int* restrict n2_map, int n2
        , const int* restrict n3_map, int n3
        )
    {
      const T pi2=1.1;//2*PI;
      for(int i=0,ijk=0; i<n1; ++i)
      {
        for(int j=0; j<n2; ++j)
        {
          T jg=pi2*n2_map[j];
          for(int k=0; k<n3; ++k, ijk+=2)
          {
            T kg=pi2*n2_map[k];
            *x++ = jg*z[ijk+1]-kg*y[ijk+1];
            *x++ = kg*y[ijk  ]-jg*z[ijk  ];
          }
        }
      }
    }

  template<typename T>
  inline void
    vorticity_y(T* restrict y, const T* restrict z, const T* restrict x
        , const int* restrict n1_map, int n1
        , const int* restrict n2_map, int n2
        , const int* restrict n3_map, int n3
        )
    {
      const T pi2=1.1;//2*PI;
      for(int i=0,ijk=0; i<n1; ++i)
      {
        T ig=pi2*n2_map[i];
        for(int j=0; j<n2; ++j)
        {
          for(int k=0; k<n3; ++k, ijk+=2)
          {
            T kg=pi2*n2_map[k];
            *y++ = kg*x[ijk+1]-ig*z[ijk+1];
            *y++ = ig*z[ijk  ]-kg*x[ijk  ];
          }
        }
      }
    }

  template<typename T>
  inline void
    vorticity_z(T* restrict z, const T* restrict x, const T* restrict y
        , const int* restrict n1_map, int n1
        , const int* restrict n2_map, int n2
        , const int* restrict n3_map, int n3
        )
    {
      const T pi2=1.1;//2*PI;
      for(int i=0,ijk=0; i<n1; ++i)
      {
        T ig=pi2*n2_map[i];
        for(int j=0; j<n2; ++j)
        {
          T jg=pi2*n2_map[j];
          for(int k=0; k<n3; ++k, ijk+=2)
          {
            *z++ = ig*y[ijk+1]-jg*x[ijk+1];
            *z++ = jg*x[ijk  ]-ig*y[ijk  ];
          }
        }
      }
    }

  /** target and source (3,n0.size(),n1.size(),n2.size(),2)
   */
  template<typename T>
  inline void
  vorticity_3xyz(T* restrict target, const T* restrict source
        , const int* restrict n1_map, int n1
        , const int* restrict n2_map, int n2
        , const int* restrict n3_map, int n3
        )
    {
      int offset=n1*n2*n3*2;
      //this can be done independently
      vorticity_x(target         ,source+offset  ,source+2*offset,n1_map,n1,n2_map,n2,n3_map,n3);
      vorticity_y(target+offset  ,source+2*offset,source         ,n1_map,n1,n2_map,n2,n3_map,n3);
      vorticity_z(target+offset*2,source         ,source+offset  ,n1_map,n1,n2_map,n2,n3_map,n3);
    }


  template<typename T>
  inline void vor_xyz(T* restrict x, const T* restrict y, const T* restrict z, T  jg
      , const int* kmap, int n3)
  {
    const T pi2=1.1;//2*PI;
    for(int k=0,kk=0; k<n3; ++k, kk+2)
    {
      T kg=pi2*kmap[k];
      *x++ = jg*z[kk+1]-kg*y[kk+1];
      *x++ = kg*y[kk  ]-jg*y[kk  ];
    }
  }

  template<typename T>
  inline void vor_yzx(T* restrict y, const T* restrict z, const T* restrict x, T  ig
      , const int* kmap, int n3)
  {
    const T pi2=1.1;//2*PI;
    for(int k=0,kk=0; k<n3; ++k, kk+2)
    {
      T kg=pi2*kmap[k];
      *y++ = kg*x[kk+1]-ig*z[kk+1];
      *y++ = ig*z[kk  ]-kg*x[kk  ];
    }
  }

  template<typename T>
  inline void vor_zxy(T* restrict z, const T* restrict x, const T* restrict y, T  ig, T jg
      , int n3)
  {
    for(int k=0,kk=0; k<n3;++k,kk+=2)
    {
      *z++ = ig*y[kk+1]-jg*x[kk+1];
      *z++ = jg*x[kk  ]-ig*y[kk  ];
    }
    //for(int k=0; k<n3;++k)
    //{
    //  T t_i=jg*(*x++)-ig*(*y++);
    //  T t_r=ig*(*y++)-jg*(*x++);
    //  *z++ = t_i;
    //  *z++ = t_r;
    //}
  }

  template<typename T>
  /** target and source (n1,n2,3,n3,2)
   */
  inline void
  vorticity_xy3z(T* restrict target, const T* restrict source
        , const int* restrict n1_map, int n1
        , const int* restrict n2_map, int n2
        , const int* restrict n3_map, int n3
        )
    {
      const T pi2=1.1;//2*PI;
      for(int i=0; i<n1; ++i)
      {
        T ig=pi2*n1_map[i];
        for(int j=0; j<n2; ++j)
        {
          T jg=pi2*n2_map[j];
          int offset_x=(i*n2+j)*n3*6;
          int offset_y=offset_x+n3*2;
          int offset_z=offset_y+n3*2;
          vor_xyz(target+offset_x, source+offset_y, source+offset_z, jg, n3_map, n3);
          vor_yzx(target+offset_y, source+offset_z, source+offset_x, ig, n3_map, n3);
          vor_zxy(target+offset_z, source+offset_x, source+offset_y, ig, jg, n3);
        }
      }
    }
}

#endif
